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ABSTRACT 

We perform two-dimensional MHD simulations of the flux emergence from 
the solar convection zone to the corona. The flux sheet is initially located mod- 
erately deep in the adiabatically stratified convection zone (—20, 000 km) and 
is perturbed to trigger the Parker instability. The flux rises through the solar 
interior due to the magnetic buoyancy, but suffers a gradual deceleration and 
a flattening in the middle of the way to the surface since the plasma piled on 
the emerging loop cannot pass through the convectively stable photosphere. As 
the magnetic pressure gradient enhances, the flux becomes locally unstable to the 
Parker instability so that the further evolution to the corona occur s. The second- 



step n onlinear emergence is well described by the expansion law by 



Shibata et al. 



( 1l989l ). To investigate the condition for this 'two-step emergence' model, we vary 
the initial field strength and the total flux. When the initial field is too strong, the 
flux exhibits the emergence to the corona without a deceleration at the surface 
and reveals an unrealistically strong flux density at each footpoint of the coronal 
loop, while the flux fragments within the convection zone or cannot pass through 
the surface when the initial field is too weak. The condition for the 'two-step 
emergence' is found to be 10^^ - 10^^ Mx with 10^ G at ^ = -20, 000 km. We 
have some discussions in connection with the recent observations and the results 
of the thin-flux-tube model. 



Subject headings: methods: numerical - MHD - Sun: chromosphere - Sun: corona - 
Sun: interior - Sun: photosphere 
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Introduction 



Solar active regions are generally thought to be the consequence of the flux emergence, 
i.e., dynamo-generated m agnetic fluxes risen from the deep convection zone due to the 



magnetic buoyancy force ((Parker 



19551 ). Observations indicate that the flux tube should 
be in a coherent form and strong enough so as not to be disintegrated by the turbulent 
motion during its emergence through the convection zone. For generating such a strong 



flux, the sub-adiabatically stratifled overshoo t region at t 



has been suggested to be the suitable place (jParker 



le bottom of the convection zone 



19751 ). Therefore, a flux emergence 



should be understood as a whole process from the base of the convective layer to the upper 
atmosphere through the surface. 

In last two decades, various series of numerical experiments have been carried out to 



investigate the physics of the flux^ 
the pioneering work was done by 



emergence. 



Shibata et al. 



or the local evolution above the surface, 



(119891 ). who performed two-dimensional 



(2D) magnetohydrodynamic (MHD) simulations of flux emergence through the undular 
mode of magnetic buoyancy instability {k \\ B, where k and B denote the wavenumber 
and the initial magnetic fleld vector, respectively; the Parker instability) to reproduce some 
dynamical features such as the rise motion of an arch fllament system and downflows along 
the magnetic field lines. Since then, the evolution of an emerging fiux and t he interaction 



with pre-existing coronal fields have been studied by 2D and 30 simulations (INozawa et al. 



1992 



20081). 



Matsumoto et al. 



Nozawa et al. 



1993 



Yokoyama fc Shibata 



1995 



1996 



Fan 



2001b 



Cheung et al. 



( 119921 ) performed the emergence from the convectively unst a 



Yokoyama fc Shibatal (11995 



3le 



19961) 



solar interior (the convective-Parker instability), while 
studied the reconnection between the expanding loop and the pre-existing fields in the 
corona and the subsequ ent formation of X-ray jets. Three-dimensional calculations by 



Matsumoto et al. 



(119931 ) were performed for the studies of the interchange {k B) and 
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quasi-interchange [k \\ B with fcifph -C 1, where ifph is the photospheric pressure scale 



height) mode instabihties. 



FanI (j2001b|) co mpared her 30 simu lation results with observed 



features of newly emerged active regions. 



Cheung et al. 



(120081 ) found that the numerical 



modeling of emerging flux regions by 3D radiative MHD simulations exhibits photospheric 
characteristics that are comparable with the observations from the Hinode/ Solai Optical 
Telescope (SOT). These simulations assumed that the initial flux is embedded just below 
the photosphere (> —2000 km) a priori as a horizontal sheet or a twisted tube. It is because 
their experiments mainly aimed to clarify the local behaviors in the solar atmosphere. Such 
an initial structure, however, is not obvious since the initial form depends on its history 
of pushing through the convection zone. Therefore, a numerical experiment including the 
evolution from a substantial depth has been needed. 

The simulations focusing on the flux emergence within the relat ively deep so lar 



interior have been done by using the thin-fiux -tube approximation (jSpruit 



D'Silva &: Choudhuri 



1993 



Caligari et al 



1981 



19951 ). One of the most important conclusions 



obtained from the various thin-fiux-tube simulations is that rising tubes with small magnetic 
flux (below 10^^ Mx for 10^ G at the base) cannot reach the photo sphere because the 



apice s of the loops loose magnetic fields and subse quently 



9951) . By assuming the anelastic approximation (iGough 



'exp 



1969 



ode' (Moreno 



Lantz fc Fan 



nsertis et al. 



199% 



Fan 



(j2001a[ ) computed the evolution of the 3D undulatory instability of the horizontal magnetic 
flux layer, which formed arch-like magnetic tubes with downflows from their apices to the 
troughs. Note that both types of approximations (thin-fiux-tube and anelastic) are not 
applicable in the upper convection zone (> —30 Mm) where the diameter of the flux tube 
exceeds the local pressure scale height and the flow velocity becomes close to the sound 
speed. 



Abbett fc Fisherl (120031 ) calculated a flux emergence by connecting the anelastic MHD 
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convective layer and the fully compressible MHD solar atmosphere from the photosphere to 
the low corona. However, their full MHD atmosphere did not include the upper convection 
zone. To investigate the detailed behavior of the emerging flux from the convection zone to 
the atmosphere, we have to deal with the full MHD numerical box including the convective 
layer, the photosphere/chromosphere, and the corona. 

The dynamical behavior of an emerging flux including both the interior and the upper 



atmosphere is not yet clear, but is thought to obey the 



'two-step emergence' model (cf. iMatsumoto et al. 



1993 



bllowing pictu re, which we name 



Magara 



20011). Magnetic fluxes 



emerged from the bottom of the convection zone are depressed and decelerated by the sub- 
adiabatic photosphere and extended horizontally around the photosphere/chromosphere. 
Meanwhile, fluxes are still transported from below to enhance the magnetic pressure 
gradient. Finally, the fluxes above the photosphere become unstable to the magnetic 
buoyancy instability so that the further evolution into the corona occurs. By the 
development of supercomputers, full MHD simulations in the convection zone without 
approximations have come to be realizable. The purpose of this work is to investigate the 
'two-step emergence' model numerically. 



I 

Several experiments have confirmed this 'two-step' model. iMagaral (120011 ) studied the 



emergence of the magnetic flux tube from the convection zone by means of 2.5-dimensional 
MHD simulations focused on the cross section of the tube. He found the deceleration of the 



rising flu x tube due to t 
outflow 



leconvectively stable photosphere and the subsequent horizontal 



Archontis et al. 



(120041 ) performed 3D simulations using the criterion by 



Acheson 



(119791) to analyze the mag netic buoyancy instability within the photosphere/chromosphere. 



while 



Murray et al. 



(120061 ) did parameter studies of the dependence of the initial magnetic 
field strength of the tube and its twist, finding that the tube evolves in the self-similar 
way when varying the field strength and that the magnetic buoyancy instability and the 
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second-step evolution do not occur when the field is too weak. In these studies, the initial 
fiux tubes are located in the uppermost convective region at some 1000 km depth. Rec ently, 



Solar Optical Telescope (SOT) aboard the Hinode satellite (e.g. 



Tsuneta et al. 



20081) 



observed small-scale magnetic fiux emergence. Tha nks to high-resolution and high-cadence 



multi-wavelength observations. 



Otsuji et al 



(I2OIOI ) found the deceleration of the apex of 



the arch filament system in the chromosphere, which can be the possible evidence of the 
two-step model. 

In this study, we perform two-dimensional fully compressible MHD simulations to 
investigate the dynamical evolution from the adiabatically stratified convective layer 
into the corona through the isothermal (strongly sub-adiabatic, i.e., convectively stable) 
photosphere/chromosphere. We set the initial magnetic fiux sheet in the moderately 
deep convection zone at 20 Mm depth, not at the bottom of the solar interior, because 
the emergence from the base is beyond the computation ability. The numerical results 
reproduce the two-step model well. Ho wever, the pic t ure of the two-step e r aergence is 



much far from the previous studies (e.g. 



Magara 



2001 



Archontis et al. 



2004 



Murrav et al. 



20061 ). The location of the initial fiux is so deep that the typical wavelengths of the Parker 



instability (10 — 20 times the local pressure scale height) are different between the primary 
and the secondary emergence; the first-step evolution occurs at z = —20 Mm and its typical 
wavelength is about 100 Mm, while, at the photosphere, the second-step emergence has 
its wavelength of the order of a few 1000 km. We also discuss the dependence of the fiux 
sheet's behavior ('two-step emergence', 'direct emergence' or 'failed emergence') on its 
magnetic field strength and amount of fiuxes through the parameter survey. 

The numerical setup and the assumed conditions used in this study are presented in 
Section O We show the typical case of the 'two-step emergence' in Section [3l while Section 
m gives the results of the parameter survey and make some discussions with preceding 
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studies and recent observations. Finally, in Section [5l we summarize the study. 



2. Numerical Model 

2.1. Assumptions and Basic Equations 

We consider an isolated magnetic flux sheet located in a convectively marginally 
stable gas layer in a two-dimensional Cartesian coordinate system (x, z), where ^-direction 
is antiparallel to the gravitational acceleration. We solve adiabatic two-dimensional 
{d/dy = 0, By = 0, Vy = 0) magnetohydrodynamic (MHD) equations. The basic equations 
are as follows; 

^ + V.(pV)=0, (1) 



BB B' 



pg = 



(2) 



OB 



(3) 



and 



{pU + p + \pV^)V + X B 
2 47r 



P9-V = Q 



(4) 



U 



1 p 



7 - Ip 



(5) 



E 



— V X B 



(6) 



where U is the internal energy per unit mass, I is the unit tensor, g = (0, 0, —go) {go is 
constant and its value is given below) is the gravitational acceleration, and other symbols 
have their usual meanings. We assume a ratio of specific heats, 7, of 5/3. 
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2.2. Initial Conditions 

The initial magnetostatic gas layer is composed of three regions: hot and cold isothermal 
layers in upper and middle regions represent a very simplified model of the solar corona 
and photosphere/chromosphere, and a non- isothermal layer in the lower region models the 



convection zone (see Fig. 1(a)). We take 2 = to be the base height of the photosphere. 
The units of length, velocity, time, and density in the simulations are Hq, Cso, Hq/Cso = To, 
and po, respectively, where Hq = k^TQlimgo) is the pressure scale height {k^ is Boltzmann 
constant and m is mean molecular mass), Cso the sound speed, and po the density at 
z = in the photosphere. The gas pressure, temperature, and magnetic field strength are 
normalized by the combinations of the units above, i.e., po = PoC'sc ^0 = '"^C'so/(7^b)5 and 
Bq = (poC^qY^^, respectively. The gravity is given as go = C^o/(7-^o) by definition. For the 
comparison of numerical results with observations, we use Hq = 200 km, Cso = 8 kms~^. 
To = Ho/Cso = 25 s, and po = 1.4 x 10^^ g cm^'^, which are typical values for the solar 
photosphere and chromosphere. Then, po = 9.0 x 10^ dyn cm~^, Tq = 4000 K, and 
Bo = 300 G. 

The initial temperature distribution of the photosphere/chromosphere and the corona 
{z > 0) is assumed to be 

T{z) = Tph + (Tcor - Tph) {tanh [{z - Zcor)/wtr] + 1} /2 , (7) 

where Tcor and Tph are the respective temperatures in the corona and in the photo- 
sphere/chromosphere, Zcor is the height of the base of the corona, and Wtr is the temperature 
scale height of the transition region. We take Tcor = IOOTq, Tph = Tq, Zcor = WHo, and 
Wtr = 0.5Ho. As for the initial temperature distribution in the convective layer (-2 < 0), we 
assume 

T{z) = Tph - az 



dT 



dz 



ad 
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where 



dT 
dz 



7-1 mgQ 



(9) 



ad 



7 

is the adiabatic temperature gradient, and a is a parameter of the temperature gradient of 
the convection zone. In our calculations, a is set to be unity, i.e., the initial temperature 
distribution in the convection zone is adiabatic. 

The magnetic field is initially horizontal, B = 0, 0), and is embedded in the 

convection zone. The distribution of magnetic field strength is given by 



B,{z) = [%Mz)/P{z)] 



1/2 



(10) 



where 



and 



P{z) - ^/f{z) , 



tanh ( 1 + 1 



- tanh ( - — — 1 + 1 

Wi 



(11) 



(12) 



Here ^5* is the ratio of gas pressure to magnetic pressure at the center of the magnetic flux 
sheet, and Zq and Zi = Zq + D are the heights of the lower and upper boundaries of the flux 
sheet, where D is the vertical thickness of the sheet. We use zq — — 100i?o — —20, 000 km. 
In all of our calculations, wq and wi are set to be 0.5i?o- We take ^* = 160 and 
D = 5Ho ~ 1000km for case 1 (the typical model), so as the initial magnetic held strength 
Bx to be 10"^ G and the total magnetic flux $ to be 10^^ Mx. To calculate the total magnetic 
flux, we regard the initial flux sheet as a rectangular prism with a base D (in z direction) 
by lOD (in y direction). 

On the basis of the initial plasma (3 distribution mentioned above, the initial density 
and pressure distribution are calculated numerically by using the equation of static pressure 
balance: 



d 
dz 



p + 



Stt 



(13) 
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The initial temperature, density, pressure, and magnetic field strength distributions for the 



typical case are shown in Fig. 1(b) 



In order to trigger the Parker instability (jParker 
density perturbations of the form 

5p = Af{z)p{x, z) cos (27rx/ A) 



19661 ) in the convection zone, small 



(14) 



are initially reduced from the magnetic flux sheet {zq < z < zi) within the finite horizontal 
domain (— 3A/4 < x < 3A/4), where A(= 400//o) is the perturbation wavelength, and 
A{= 0.01) is the maximum value of the initial density reduction. The definition of f{z) is 
given in equation f|T2|) . 



2.3. Boundary Conditions and Numerical Procedures 

The domain of the simulation box is {x^ain < x < Xmax) and {z^^^ < z < -Zmax), where 
Xrain = — 400ifo! x-^aa.^ = 400ifo) ^min = — 200ifo) ^Jid z^^-^ = 200Hq, i.e., the total size of the 
box is 160 Mm x 80 Mm. This is much larger than those of the calcul ations focusing on 



Shibata et al. 



19891 . 



the emergence from the uppermost convection zone to the corona (e.g. 
etc.). Periodic boundaries are assumed for x = Xmm and x = Xmax, symmetric boundaries 
for z = z^i^ and z = Zmax- A damping zone is attached near the top boundary to reduce 
the effects of reflected waves. 

To solve the equations numerically, we use the modified Lax-Wendroff scheme. The 
code is a part of the numerical package CANS (Coordinated Astronomical Numerical 
Software) maintained by Yokoyama et alQ . For the typical model (case 1), the total number 



^ CANS (Coordinated Astronomical Numerical Software) is available online at: 



http : //www- space . eps . s .u-tokyo . ac . jp/$\sim$yokoyama/etc/cans/ 
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of grid points is (A^^; x A^^) = (1536 x 1920), and the mesh sizes are Ax = 0.52Ho and 
Az = 0.21ifo) both of which are uniform. For comparison, we perform other simulations 
with different values of the parameters. We vary the values of the initial magnetic field 
strength and of the total magnetic flux $ of the emerging flux sheet by adjusting 13^, and 
D. In these calculations we set the total number of grid points [N^ x A^^) = (1024 x 1280), 
and the uniform mesh sizes Ax = 0.78iJo and Az = O.SIHq. The cases we examine are 
summarized in Table [TJ 



3. General Evolution 

In this section we show the numerical results of the typical case (case 1: = 10^ G 
and $ = 10^^ Mx); the results display the 'two-step emergence'. Figure [2] illustrates the 
development of the density profile with magnetic field lines and velocity vectors, while 
Figure |3] indicates the height of the apex of the magnetic field Zapex and its rise velocity 
Vzapex- The time evolution can be divided into four stages according to the rise velocity at 
the apex of the emerging field. In the first stage (0 < t/ro < 700), the magnetic flux begins 
an emergence in the convection zone due to the Parker instability driven by the magnetic 



buoyancy (Fig. 2(b)). The rise velocity increases continuously in this stage (Fig. [3]). The 
second stage (700 < t/ro < 1900) is characterized by a gradual deceleration (Fig. |3]) because 
the arch-like emerging field becomes deformed to horizontal so that the mass on the apex 
can no longer fall down along the field lines, and thus continues to pile up on the horizontal 



field. As a consequence, the emerging flux stretches around the solar surface (Figs. 2(c 



and (d)). In the third stage (1900 < t/ro < 2000), the top part of the emerging magnetic 
flux almost stops at the surface while fluxes are still emerging from below, that is, the 
magnetic pressure gradient on the upper boundary of the flux sheet continues to enhance. 
When the magnetic pressure gradient gets steepened enough, the Parker instability sets in 
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and drives the further evolution into the upper atmosphere (Figs. 2(e) and (f)). In the 



final, fourth stage {t/ro > 2000), the magnetic fiux evolves to the corona d ue to magneti c 



press ure, which is consistent with the results of classical calculations (e.g. 



Shibata et al. 



19891 ) (Figs. 2(g) and (h)). In the following, we will discuss each stage in more detail, and 
examine the dynamical structure of the expanding loop and the related forces acting on the 
magnetic field. 



3.1. First Stage (0 < t/ro < 700) 

The initial sinusoidal perturbation in the fiux sheet triggers the Parker instability so 
that the fiux sheet begins to rise through the solar convection zone by magnetic buoyancy 



(Fig. 2(b)). In this phase, the rise speed enhances continuously as the fiux sheet emerges. 



The rising fiux becomes arch-like shape owing to the stronger buoyancy of the loop center. 
The evacuation in the apex by the downfiow along the field lines due to gravity leads to 
the acceleration of the loop until reaching the local Alfven speed (Ca = B/^/A^Tp). Figure 
m is a close-up view of the evolution between t/ro = 400 and t/ro = 600 of Figure [3l 
where a dashed line indicates the rise velocity of the apex, while the local Alfven speed is 
represented by a solid line. This figure shows that the rise velocity increases so as to get 
close to the local Alfven speed. 



3.2. Second Stage (700 < t/ro < 1900) 



At around t/ro = 700, the rise velocity changes from acceleration to deceleration (Fig. 



[3]), and at t/ro = 1000, both legs of the loop become vertical (Fig. 2(c)). The central part 



of the emerging loop fiattens and expands horizontally along the surface. Drained gas from 
the apex fiows down along the field lines to each trough, so that the both troughs sink into 
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the deep convection zone. 



In this stage, the rise motion turns into the deceleration as seen from Figure El iMagara 



(I2OOII ) found that the tube reduces its rise speed and becomes flatt ened since the rise 



Murray et al. 



fcooeh 



motion cannot persist through the convectively stable photosphere, 
explained that the deceleration process occurs because the downward pressure gradient 
exceeds the upward magneti c buoyancy w h en th e emerging flux tube is close enough 



to the surface. In addition, 



Murray et al. 



( 120061 ) found a period when the rise speed 



diminishes due to the aerodynamic drag exerted by the flows surrounding the tube while 
in the convection zone. In our case, however, the deceleration occurs at much deeper level 
{z ~ —50Hq = —10,000 km) than the previous studies {z ~ —850 km). Moreover, our 
simulation is carried out in a two-dimensional scheme, while the previous studies were done 
in 2.5D or 3D, so that a three-dimensional force such as the aerodynamic drag does not 
exert on our emerging loop. Our results are explained by another mechanism as follows. 

Figure |5] shows the distribution of the density difference from the initial state 
(Ap = p{t) — p(0)) and the horizontal component of the magnetic field (5^) along the z-axis 
at t/TQ = 600, t/ro = 1000, and t/rg = 1960. As can be seen in Fig. [5l the mass piled on the 
emerging loop cannot rise through the photosphere/chromosphere ranging from z/Hq = 
to z/Hq = 10 because this isothermal (i.e. strongly sub-adiabatic) layer is convectively 



stable. Figure 6(a) shows the vertical component of the forces acting upon the apex of 
the loop. They are gas pressure gradient —Vp = —dp/dz, gravity pg^, magnetic pressure 
gradient — Vpmag = —d[B^/ {S7T)]/dz, and magnetic tension tma.g = [{B ■ V)-B]^/(47r). Figure 



6(b) shows the acceleration calculated by dividing the total force with the gas density 
{Fz/p = {—Vp — pgo — Vpmag — tmag)/p)- It should be noted that the total force is much 
smaller than each force since the rising loop is almost in a mechanical equilibrium. Figure 



6(b) indicates that the acceleration turns from positive to negative at around t/rQ = 870, 
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which means the rise velocity changes into deceleration phase at that time. 



The deceleration of the crest and the continuous rise motion of the both sides cause 
the loop flattened, which, in turn, makes the mass left on the flattened loop. As a result, 
the rising flux decelerates and stretches horizontally beneath and around the surface. This 



process can possibly exp^ 
much smaller scales (e.g. 



ain the formation o f the 'initial flux' of the previous studies in 



Shibata et al. 



1989 



etc.). We have to carry out three dimensional 
experiments because another effects such as the aerodynamic drag would act on the actual 
expanding loops. However, the above-mentioned scheme can explain the deceleration in the 



convection zone if the emerging field has a sheet-like shape rather than a tubular form. 



3.3. Third Stage (1900 < t/ro < 2000) 



The field is decelerated and flattened due to the isothermal (sub-adiabatic) stratification 
on the surface; meanwhile, the fluxes are continuously transported from beneath, that is. 



the magnetic pressure gradient keeps enhancing at the surface (Fig. 2(e)). At the location 
where the magnetic pressure gradient gets steepened enough, the further evolution to 
the corona occurs on the ba sis of the Parker instability. At the point of the second-step 



emergence, a 'pressure hill' (lArchontis et al. 



2004; 



Magara 



200 ll ) is formed, which indicates 



that the plasma in the photosphere drains out sideways and the magnetic field covers this 
area, i.e., the stratification is top-heavy (see Fig. [7]). 

To confirm the onset of the Parker instability, we use the criterion ob tained by 



NewcombI (jl96ll ). The criterion for the Parker instability is (jNewcomblll96ll ) 

dp 



dz 7p 



We plot the index 



, dp p^Qq 

ijj = , 

dz "jp 



(15) 



(16) 
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that is, the area with negative ip is subject to the instabihty. Figure [8] illustrates the ip 
distribution with field lines just before the second-step emergence at t/ro = 1960. This 
figure indicates that the index is negative at around the point of emergence; we can 
conclude that the further evolution to the corona is ascribed by the Parker instability. 

At t/TQ = 1960, plasma /3 = p/{B^ /Sn) is order of unity (~ 2) and the magnetic 
field strength is about 700 G at t he low photosphere wi thin the emergent area, which are 



consistent with observations (e.g. 



Watanabe et al 



20081) 



3.4. Fourth Stage (t/ro > 2000) 



In this final phase, the magnetic fiux emerging within the photosphere begins to expand 
to the solar corona by the magnetic pressure on the condition that the gas pressure acting 



on the surface of the fiux is weak enough (Figs. 2(f) (h)). The expanding loop finally forms 
a lar ge coronal loop. Th is process is similar to that of the results of classical calculations 



[e.g. 



Shibata et al. 



3 



19891, etc.). 



The characteristics of this nonlinear phase is a self-similar evolution. Figures 9 (a 



indicate the distributions of the vertical component of velocity, gas density, and horizontal 



magnetic field strength along the axis x/Hq = 20, where the apex of the loop is 



Figs. [2(g)] and [(h)]), at t/ro = 2000, 2020, and 2040. According to 
expansion law is written as 



Shibata et al. 



ocated (see 



(119891), the 



sO 



p (X Z 

B,. oc z-^ 



az/Ho, 

4 



(17a) 
(17b) 
(17c) 



where a = 0.05 — 0.07 is about half the non-dimensional linear growth rate for plasma 
P = 0.5 — 2.0 of the fiux sheet. In our study, plasma /3 of the magnetic field has been 
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calculated to be 2 at I/tq = 1960 before the further evolution begins (see Section [373]) . 
which suggests that the velocity- height relation is 

VjC,or^0.05z/Ho. (18) 



This relation is overplotted in Figure 9(a) with a solid line. The other relations given by 



equations (]17bp and fll7cl) are also overplotted in Figures 9(b) and (c) with solid lines. As 
seen from Figure [9l the nonlinear growth to the solar corona is consistent with that of 



Shibata et al. 



(119891 ). 



The size of the coronal loop at t/ro = 2070 is found to be of AOOHq = 80, 000 km width 
and 200/^0 = 40, 000 km height, while, at the surface, plasma /3 is of the order of unity and 
the field strength is about 4Eo(= 1200 G). 



4. Parameter Survey and Discussion 

We carry out a parameter survey by changing the values of the initial field strength 
Br^ and the total magnetic flux $ of the emerging flux sheet to investigate the condition of 
the sheet's behavior. A summary of the values of and $ under consideration is given in 
Table [H Figure [10] shows the results of the parameter survey, where diamonds, asterisks, 
and X's represent 'two-step,' 'direct,' and 'failed' evolution, respectively. Fluxes belong 
to the direct emergence group do evolve into the corona, however, they do not reveal the 
deceleration, unlike fluxes of the two-step emergence, while those of the failed emergence 
fragment within the convection zone or cannot pass through the photosphere. Figure [TT] 
indicates height-time relations of the fluxes along the axis x/Hq = 0. In this section, we 
will discuss each group in detail. 
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4.1. Direct emergence 

In cases 2 and 3, the fluxes show the emergence to the corona without any deceleration 
at the surface. In other words, they 'directly' emerge to the corona. As shown in Figure 



11(a) , the height-time curves of cases 2 and 3 do not have an inflection point. We name them 
the 'direct emergence' group. The absence of an inflection implies that their evolutions are 
not affected by the isothermal photosphere/chromosphere at all since they have extremely 
strong field (note that these cases can be found in the upper right of Fig. [TU]) . 

Figure [12] shows the time evolution of the density profile, magnetic field lines, and 
velocity vectors for case 2. Color contour is the same as that of Fig. [2l Each flux of 
this group exhibits field strength B ~ 10^ G and plasma /3 ~ 0.1 at the surface after the 
emergence, which is not consistent with observations. Therefore, 'direct emergence' model 
is not suited for the formation model of active regions. 



4.2. Two-step emergence 

Cases 4-7 show the 'two-step emergence' to the corona as well 1 (typical model). 



Height-time relations of this group are shown in Fig. 11(b) Each line of this group has 
an inflection point beneath the surface {z/Hq = 0), that is, rise velocity of the emerging 
flux turns from acceleration to deceleration phase within the convection zone due to the 
isothermal photosphere/chromosphere. As the figure indicates, the temporal length of the 
rising stage decreases with increasing initial field strength and total flux. 

Among these five 4, 5, and 6 show unrealistically strong flux densities 

B ~ 10^ G (plasma (3 ~ 0.1) at the photosphere after the coronal loops are built up. The 
realistic models of the formation of active regions are cases 1 and 7 ($ = 10^^ — 10^^ Mx 
with Bx = 10^ G at z = —20,000 km). This result gives an important suggestion that the 
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fluxes which form active regions are hkely to have experienced the 'two-step emergence. 



Our numerical results agree with recent satellite data. lOtsuji et al.l ( 120101 ) found the 
deceleration and the horizontal spreading of an emerging flux within the solar chromosphere 
by using Hinode/ SOT. This observation supports the concept of the 'two-step emergence' 
model. At the same time, there is a difference between our results and the observation. 
The deceleration occurs in the chromosphere, not beneath the surface. The difference may 
partly come from the structure of the emerging loop. The numerical flux considered here 
has a sheet-like structure, on which the plasma piles continuously during its evolution (see 
Section [3]). If the emerging loop is part of a (twisted) flux tube, the plasma on the loop 
can drain around the cross-section of the tube. Therefore, the emerging flux tube may 
rise through the convection zone faster than the sheet-like flux, and the tube may suffer a 
deceleration at a higher alt itude when the tube itself passes through the convectively stable 



layers. On the other hand, iMagaral (l200ll ) reported that even in case with an initial twisted 
tube, there is a deceleration close to the photosphere. The attribution of the deceleration 
altitude to the structure of the magnetic flux is still oversimplified at this moment. More 
quantitative study by three-dimensional simulations concerning with this issue is necessary. 

It is possible that the emerging flux suffers a deceleration beneath the surface and 
exhibits the 'two-step' evolution, or, in some cases, fluxes show the 'multi-step' emergence 
since the structure of the rising loop in the convection zone cannot be achieved by the 
optical observations. To know the structure and the behavior of the flux emerge nce below 



the photosphere, advanced local helioseismology is needed (e.g. 



Sekii et al 



2003) 
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4.3. Failed emergence 

Cases given in Fig. [10] with X's belong to the 'failed emergence' group (cases 8-14). 
Fluxes of cases 9-14 suffer a fragmentation due to the continuous motion within the 
convection zone, so that further emergence does not occur (see Fig. [T^ . It is because the 
fluxes have weak fields compared to the local kinetic energy density of flow motion induced 
by the initial perturbation. Figure [T3] shows the ratio of the magnetic energy density 
-E'mag to the local kinetic energy density i^kin of case 14 along x/Hq = at t/rp = 4000 



(Fig. 13(e)), where -Emag = B /(Sir) and _Ekin = pv /2, respectively. It reveals that the 



magnetic energy is weaker than the kinetic energy all over the convective layer. Height-time 



relations of these fluxes are presented in Figure 11(c) As the figure shows, each line of this 
group reveals a continuous fluctuation and remains around the surface, indicating that the 
continuous motion repeats within the solar interior. On the other hand, flux of case 8 do 
not show the secondary emergence because the flux fails to enhance the magnetic pressure 
gradient while crossing the surface. This flux keeps its coherency all the while since the 



field strength is almost in the same range as the local kinetic energy density (see Fig. [15] 
Dotted line in Fig. [TU] indicates the criteria f or the 'explosion' of the eme rging flux 



obtained by thin- flux-tube (TFT) approximation (iMoreno-Insertis et al.l 119951 ). They 
found that the rising tubes with a small amount of flux cannot reach the surface due to 
the 'explosion': if the tube rises sufficiently slowly, the stratification inside the tube gets 
close to a hydrostatic equilibrium along the field lines while the stratification outside is 
super-adiabatic. When the pressure difference between inside and outside the tube is small 
enough at the base of the convection zone, the difference decreases as the tube rises because 
the pressure gradient inside the tube is less steeper than that outside, so that the magnetic 
field at the apex can no longer be confined at a certain 'explosion' depth. According to their 
calculations, the tubes with the initial field strength and magnetic flux (10^ G, 10^^ Mx) 
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and (10^ G, 10^^ Mx) at the base of the convection zone can reach close to the surface. 
Their field strength at t/Rq = 0.97 {z = —20, 000 km, i.e., the depth of ou r initial fluxes) 



Moreno- Insertis et al 



1995 



Fig. 1). We 



are 3 x 10^ G and 2 x 10"^ G, respectively (see 
adopt them as criteria for the 'non-explosion'. The majority of the fluxes with parameters 
in the range where they could not have reached at z = —20, 000 km according to the TFT 
model (left to the dotted line) also ceases emergence even in our MHD model. 

There are, however, some cases in which, although the fluxes are expected to reach at 
z = —20,000 km level by the TFT model, they fail to evolve further (cases 8, 9 and 10). 
Especially, case 8 (10^ G with 10^° Mx at 2; = —20, 000 km) reveals the photospheric field 
strength _B ~ 1 G, indicating it could be the source of the magnetic fields in the quiet sun 
if the fields are enhanced, e.g., by the flux expulsion due to the magneto-convection at the 
surface. Interestingly, this flux maintains its coherency all the while and remains floating 
around beneath the surface after it fails the further evolution to the corona because it has a 
strong field that is in equipartition with the kinetic energy density (Fig. [T5|) . This suggests 
that the flux of case 8 can be the origin of the epheme ral regions. Flu x 8 has 10^" Mx, 



consistent with the flux of mid-sized ephemeral regions (IHagenaar 



20011). 



Conclusions 



We perform the nonlinear two-dimensional simulations to investigate the behavior of 
emerging flux from moderately deep convection zone {z = —20,000 km). We set a much 



wider numerical box ( 
Parker instability (e.g. 



60 Mm X 80 Mm) than those of the previous experiments on the 



Shibata et al. 



19891 ) 



In the typical case {B = 10^ G with $ = 10^^ Mx at 2 = -20,000 km), the results 
show the 'two-step emergence'. In the middle of the way of the first emergence to the solar 
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surface, the flux loop turns from an acceleration phase to deceleration when approaching 
the (sub-adiabatically stratified, i.e., convectively stable) photosphere/chromosphere. The 
emerging flux has a sheet-like shape, thus it is difficult for the mass on the loop to escape 
from the area between the loop and the convectively stable surface. This mass pile-up 
causes the loop decelerate. The deceleration of the apex of the expanding flux and the 
continuous rising of the hillsides make loop flattened, which results in the plasma kept on 
the flux. This deceleration mechanism is another new one and different from those of the 



preceding studies (IMagara 



2001 



Archontis et al. 



2004 



Murray et al 



20061 ) with magnetic 



flux tubes in much smaller regions. However, our result predicts the behavior of a flux 
within the convection zone, provided the flux has a sheet-like structure . As a result of 
the deceleration and the flattening, the flux spreads sideways just beneath the surface, at 
which point the rise velocity of the crest of the loop is almost zero. Meanwhile, the flux is 
continuously transported from below, then the magnetic pressure gradient enhances locally 
in the photosphere. We found that the further evolution to the corona occurs on the basis 
of the Parker instability. At the point of the instability, plasma /3 is calculated to be order 
of unity (~ 2) and the magnetic field strength is about 700 G. In the final stage, the fiux 
shows the nonlinear evolution to the corona, which resembles the classical experiments (e.g. 



Shibata et al. 



law by 



19891. e t c.). 



Shibata et al. 



'he second-step evolution is described clearly by the expansion 



(119891 ). We find that the coronal loop exhibits 80, 000 km width with 



40, 000 km height, while the field strength of each footpoint at the surface is about 1200 G. 

We perform parameter runs by changing the initial field strength S^. and the total 
fiux $ to investigate the condition of the 'two-step emergence'. The results of the runs 
under considerations can be divided into three groups: 'direct', 'two-step', and 'failed' 
emergence. In case of the 'direct emergence', the fiux do evolve to the corona, but they do 
not show the deceleration by the isothermal surface due to their strong initial magnetic 
fields (10^^ — 10^^ Mx with 10^ G at ^ = —20, 000 km). The coronal loops present irregularly 
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strong flux densities at the footpoints; thus, we conclude that they are not suitable for the 
formation models of active regions. As for the cases showing the 'two-step emergence', two 
out of five exhibit the favorable values of the photospheric field strength and plasma /3. The 
others have so large values that they cannot be regarded as realistic models of active regions. 
We can say that active regions on the sun are likely to have undergone the deceleration and 
likely to show the 'two-step emergence' mentioned above. The condition for this 'two-step' 
active region is ranging from 10^^ to 10^^ Mx with 10^ G at z = —20, 000 km in the 
convection zone. Some recent observations support this two-step model. The cases with 
B < 10"^ G reveal 'failed' evolutions; they fragment within the convection zone or cannot 
have sufficient magnetic pressure gradient to trigger the instability that the second-step 
emergence do not occur although the fiux maintains its coherency. We have some discussions 



i n con nection with the results of the thin-fiux-tube (TFT) model by 



Moreno-Insertis et al. 



( 119951 ). The cases which are found to have 'exploded' in the deeper point in the TFT 
scheme also do not show further evolutions in our MHD scheme. However, there are some 
cases which escape the 'explosion' fail the second-step evolutions, one of which is possibly 
the source of the magnetic field in the quiet sun. 

The present calculations are in a two-dimensional scheme solving simplified equations. 
Thus we have to demonstrate the more realistic experiments in 3D. At the same time, 
advanced observations by helioseismological technique are needed to reveal the detail of the 
emerging fiux in the convection zone. 



Numerical computations were in part carried out on NEC SX-9 at Center for 
Computational Astrophysics, CfCA, of National Astronomical Observatory of Japan. 
Numerical computations were in part carried out on Space Science Simulator (NEC SX-6) 
of JAXA Supercomputer System. Numerical computations were in part carried out on 
NEC SX-8 at Nobeyama Solar Radio Observatory of National Astronomical Observatory of 
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Fig. 1. — (a) Schematic depiction of the initial set-up. (b) One- dimensional (2;-) distributions 
of the initial density (solid line), pressure (dotted line), temperature (dashed line), and 
magnetic field strength (dashed-dotted line). 
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Fig. 2. — Time-evolution of the 'two-step emergence' for case 1 (typical model), (a) t/ro = 0; 
(b)l t/ro = 500; \(c)\ t/ro = 1000; M)\ t/ro = 1500; 



t/ro = 1960; (f) t/r^ = 2020; (g 



t/ro = 2040; (h) t/ro = 2070. Logarithmic density profiles (log^o (p/po)) are indicated by 
color contour, while magnetic field lines and velocity vectors are overplotted with black lines 
and arrows. The white box at t/ro = I960 shows the area we analyze the Parker instability 
(see Fig. |8]). This figure is also available as an avi animation in the electronic edition of the 
Astrophysical Journal. 
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Fig. 3. — Time evolution of the apex of the flux sheet; sohd line is a height of the apex, whose 
rise velocity is indicated with dashed line. Photospheric level {z/Hq = 0) is overplotted with 
dotted line. 




Fig. 4. — The close-up view of Fig. [3l Upward velocity at the apex of the loop is plotted by 
a dashed line, while the local Alfven speed is represented by a solid line. 
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Fig. 5. — Density differences from the initial state (Ap = p{t) — p{0), solid lines) and 
horizontal field components {B^, dotted lines) along the z-axis at the center of the simulation 
box {x/Hq = 0) at the three different times. As the magnetic flux rises, mass piles up on the 
loop. However, the mass cannot persist through the isothermal photosphere (0 < z/Hq < 
10). 
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Fig. 6. — (a) Time evolution of the vertical components of the forces acting on the apex of 
the emerging loop: gas pressure gradient (— Vp = —dp/dz, solid line), gravity {pgo, dotted 
line), magnetic pressure gradient (— Vpmag = —d[B'^/{87r)]/dz, dashed line), and magnetic 
tension (tmag = [{B ■ V)-B]2/(47r), dashed-dotted line), (b) The same of the acceleration 
{Fz/p = (— Vp — pQo — Vpniag — ^mag)/p)- A dottcd line shows that the acceleration equals 
zero. 
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Fig. 7. — Vertical distributions of magnetic pressure (solid line), gas pressure (dotted line), 
and gas density (dashed line) along the axis x/Hq — 20 at t/ro = 1960, which is the central 
position of the second-step emergence. 
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Fig. 8. — Spatial distribution of the index i/j = —dp/dz — p'^Qo/i'jp)- at the location shown 



by the white box in Fig. 2(e), Color contour indicates ip (blue represents negative), while 
magnetic field lines are overplotted by solid lines. The aspect ratio is arranged. The index 
ip is negative around the area of the further evolution {x/Hq = 20, z/Hq = 0), indicating 
that this area is subject to the Parker instability. 
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Fig. 9. — (a) Distribution of the upward velocity along the vertical axis x/Hq = 20. Dotted, 
dashed, and dash-dotted lines indicate the distribution at t/ro = 2000, t/ro = 2020, and 
t/ rn = 2040, respective ly. A solid line shows the theoretical velocity- height relation according 



to 



Shibata et al. 



notation of lines is the same as in 



( 1l989l ). (b) Distribution of the gas density along the axis x/Hq = 20. The 

Distribution of the horizontal component of the 



magnetic field along the axis x/Hq = 20. The notation of lines is the same as in 
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Fig. 10. — Results of the parameter survey. Diamonds, asterisks, and X's represent two-step 
emergence, direct emergence, and failed emergence, respectively. Case numbers are plotted 
on the lower right of each symbol. Criteria f or the 'explosion' of the eme rging flux obtained 



by thin-fiux-tube approximation simulations (IMoreno-Insertis et al. 
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Fig. 11. — (a) Time evolution of the height of the flux sheet along the axis xjH^ = 
for 'direct emergence' (cases 2 and 3). Dashed and dashed-dotted lines represent evolutions 
of cases 2 and 3, respectively, while solid line is the height-time relation of case 1 (typical 



model), (b) Same for 'two-step emergence' (cases 1 and 4-7). Solid, dashed, dashed- 
dotted, dashed-dotted-dotted-dotted, and long-dashed lines represent evolutions of cases 1, 



4, 5, 6, and 7, respectively, (c) Same for 'failed emergence' (cases 8-14). Dashed (thick), 
dash-dotted (thick), long-dashed (thick), dashed-dotted-dotted-dotted, dashed (thin), and 
dash-dotted (thin), long-dashed (thin) lines represent evolutions of cases 8, 9, 10, 11, 12, 13, 
and 14, respectively, while solid line is the height-time relation of case 1 (typical model). In 
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Fig. 12. — Time-evolution of the 'direct emergence' for case 2. (a) I/tq = 0; (b) t/ro = 200 



t/ro = 300. Logarithmic density profiles (log^^g (p/Po)) are indicated by color contour, 



while magnetic field lines and velocity vectors are overplotted with black lines and arrows. 
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Fig. 13. — Time-evolution of the 'failed emergence' for case 14. (a) t/rg = 0; (b) t/rg = 1000 



t/ro = 2000; (d) t/ro = 3000; 



t/ro = 4000; [(f)] ^^"0 = 5000. Logarithmic density 
profiles (log^g (p/Po)) are indicated by color contour, while magnetic field lines and velocity 
vectors are overplotted with black lines and arrows. 
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Fig. 14. — The ratio of E^^^ = ^^/(STr) to £^kin = pv'^/2 of case 14 along x/Hq = at 
t/ro = 4000. 
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Fig. 15. — Time-evolution of the 'failed emergence' for case 8. (a) t/ro = 0; (b) t/ro = 1000 



t/ro = 2000; (d) t/ro = 3000; 



t/ro = 4000; [(f)] V^o = 5000. Logarithmic density 
profiles (log^g (p/Po)) are indicated by color contour, while magnetic field lines and velocity 
vectors are overplotted with black lines and arrows. 
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Table 1. Summary of Cases 



case [G]'' $ [Mx]'' /S*" D [km]'* N^, x A^/ 



1 


1 


.0 


X 


10^ 


1, 


.0 


X 


10^1 


1.6 X 10^ 


1000 


1536 


X 


1920 


2 


8, 


,1 


X 


10^ 


9, 


.8 


X 


1023 


1.1 X 10^1 


11,000 


1024 


X 


1280 


3 


1 


.0 


X 


10^ 


1, 


.0 


X 


1023 


2.0 X IQ-i 


3000 


1024 


X 


1280 


4 


1 


.1 


X 


10^ 


1, 


.0 


X 


1022 


2.0 X 10-1 


840 


1024 


X 


1280 


5 


1 


.1 


X 


10^ 


1, 


.1 


X 


1021 


2.0 X 10-1 


200 


1024 


X 


1280 


6 


1 


.1 


X 


10* 


1. 


.1 


X 


1023 


2.2 X 102 


11,400 


1024 


X 


1280 


7 


1 


.0 


X 


10* 


1. 


.0 


X 


1022 


1.6 X 102 


3200 


1024 


X 


1280 


8 


8, 


.9 


X 


103 


1, 


.0 


X 


1020 


1.5 X 102 


260 


1024 


X 


1280 


9 


1 


.1 


X 


10^ 


1, 


.0 


X 


1022 


2.2 X 10* 


11,000 


1024 


X 


1280 


10 


1 


.0 


X 


10^ 


1, 


.0 


X 


1021 


1.6 X 10* 


3200 


1024 


X 


1280 


11 


9 


.9 


X 


10^ 


1, 


.1 


X 


1020 


1.5 X 10* 


1000 


1024 


X 


1280 


12 


1 


.0 


X 


102 


1. 


.0 


X 


1021 


2.5 X 10^ 


11,400 


1024 


X 


1280 


13 


1 


.1 


X 


102 


1, 


.0 


X 


1020 


1.4 X 10^ 


3200 


1024 


X 


1280 


14 


1 


.0 


X 


102 


1, 


.0 


X 


1019 


1.4 X 10^ 


960 


1024 


X 


1280 



"^Initial field strength. 

''Total magnetic flux. 

'^Plasma beta at the sheet center. 

^Width of the sheet. 

''Total grid points. 



